Regulation of Antimicrobial Peptide Activity via Tuning Deformation Fields by Membrane-Deforming Inclusions

Antimicrobial peptides (AMPs) are considered prospective antibiotics. Some AMPs fight bacteria via cooperative formation of pores in their plasma membranes. Most AMPs at their working concentrations can induce lysis of eukaryotic cells as well. Gramicidin A (gA) is a peptide, the transmembrane dimers of which form cation-selective channels in membranes. It is highly toxic for mammalians as being majorly hydrophobic gA incorporates and induces leakage of both bacterial and eukaryotic cell membranes. Both pore-forming AMPs and gA deform the membrane. Here we suggest a possible way to reduce the working concentrations of AMPs at the expense of application of highly-selective amplifiers of AMP activity in target membranes. The amplifiers should alter the deformation fields in the membrane in a way favoring the membrane-permeabilizing states. We developed the statistical model that allows describing the effect of membrane-deforming inclusions on the equilibrium between AMP monomers and cooperative membrane-permeabilizing structures. On the example of gA monomer-dimer equilibrium, the model predicts that amphipathic peptides and short transmembrane peptides playing the role of the membrane-deforming inclusions, even in low concentration can substantially increase the lifetime and average number of gA channels.


Introduction
Antimicrobial peptides (AMP) are considered as prospective candidates for the next generation of antibiotics; they are thought to be able to overcome multiple drug resistances of bacteria. In a variety of AMPs several (overlapping) classes can be distinguished: αhelical amphipathic peptides [1,2], β-hairpin peptides [3], channel-formers [4]. AMPs can fight bacteria via several mechanisms, among which the cell lysis or pore formation in the plasma membrane is relatively widespread [1,[5][6][7]. The selectivity of AMPs towards bacteria is mainly based on negative electric charge of the outer leaflet of bacterial plasma membrane, as opposed to electrically neutral outer leaflets of membrane of eukaryotic cells. Thus, most of AMPs carry large positive electric charge to ensure strong binding of peptides to bacterial membranes [1,3,6].
Two types of pores are usually observed: (i) barrel-like, where the pore edge is formed exclusively by peptide molecules [8,9]; (ii) toroidal, where the edge is formed by both peptides and lipids [8,10,11]. The characteristic feature of amphipathic peptides having α-helical secondary structure is that their helix has one predominantly hydrophobic side surface, while the opposite side surface is mainly polar and charged. The ratio of hydrophobic/hydrophilic areas determines the so-called polar angle, which is considered to be responsible for the type of the pore formed by the peptides: larger polar angle yields barrel-like pores. AMPs having β-sheet secondary structure mostly form barrel-like pores of fixed stoichiometry, and, consequently, fixed pore diameter and conductivity [3].
The mechanisms of pore formation by membrane-bound peptides are still under debate. Experimentally, two configurations of AMPs are commonly observed: (i) a surfacebound S-configuration, in which the AMPs lie almost parallel to the membrane surface, having their hydrophobic parts incorporated into the lipid monolayer hydrophobic core; (ii) a membrane-spanning I-configuration, in which AMPs stand about perpendicularly to the membrane plane forming the edge of toroidal or barrel-like pore [8,10,12]. An intensive transformation from S to I-configuration leading to pore formation takes place at around certain so-called critical surface concentration of AMP, which usually weakly depends on particular AMP chemical structure, and for the most AMPs has the order of magnitude of about 40-100 lipids per peptide [1,8,10,13]. In the membrane-bound S-configurations AMPs incorporate their hydrophobic parts into the lipid monolayer, thus unavoidably inducing elastic deformations of the membrane [2,12,[14][15][16]. The characteristic length of lateral spread of membrane deformations is about several nanometers [17][18][19]. When the membrane-deforming peptides are far separated, their induced deformations are independent and the corresponding elastic energies are additive. However, as the distance between the peptides decreases, the deformations overlap leading to effective lateral interaction between the peptides [19]. Deformations induced by a typical α-helical amphipathic peptide affect about 80 lipid molecules in the peptide vicinity [19,20] in both membrane leaflets; this value is pretty close to the AMP critical surface concentration. In other words, an intensive pore formation starts when, on the average, the whole membrane is covered by peptide-induced deformations. Substantially higher AMP concentrations usually lead to micellation of the membrane; in this case the AMP acts as a detergent [21]. Physical reasons of cooperative re-orientation of initially membrane-bound AMP molecules in the course of transformation to I-configuration, as well as the reasons of similarity of the critical concentration values among AMPs having different chemical structures are not yet resolved. In the work [22] the authors thoroughly analyzed approximately 16,000 different peptides 9-15 amino acid residues long, soluble in water and increasing the permeability of phospholipid membranes upon binding to them. As a result, 10 peptides with high rates of solubility and membrane permeabilization were selected. The analysis of the selected peptides did not reveal any unique features of the chemical structure. These extremely active peptides had no common sequence motifs, differed in peptide length and total electrical charge. However, the peptides had a similar secondary structure and core hydrophobicity. In an extensive review on membrane permeabilizing peptides [23], it is noted that there is no established single mechanism of action of such peptides, and that the membrane permeabilization is the result of a heterogeneous and dynamic ensemble of structural states of peptides that differ depending on a variety of experimental conditions. In this regard, there are practically no discrete active structures among the thousands of known AMPs, and there are no clear rules for determining the chemical structure of a peptide to achieve a given membrane activity. In this review, the authors propose not to analyze specific features of the chemical structure of peptides, but rather to consider the distribution of the totality of all membrane-permeabilizing peptides in the space of their mechanical action on lipid membranes.
Unfortunately, at the bulk concentrations of the AMP that are necessary to achieve the critical surface concentration, most AMPs effectively bind and permeabilize membranes of eukaryotic cells as well [1]. The AMP affinity towards a membrane of particular lipid composition is determined by peptide electric charge and by the tendency to hide the hydrophobic parts of the molecule inside the membrane hydrophobic core. The latter is achieved via partial incorporation of the peptide into the lipid monolayer, and it is the incorporation that is responsible for membrane deformations leading to peptide-peptide lateral interaction and cooperative pore formation. Thus, an effective AMP has to have a substantial hydrophobic part of the molecule. The hydrophobic interactions depend but weakly on the particular type of the membrane (bacterial or eukaryotic), and for this reason there should be always a fraction of AMPs that bind to eukaryotic cells, potentially leading to their permeabilization or even lysis. Such concomitant damage of eukaryotic cells substantially hinders the development and practical use of AMPs. A way to overcome this side effect is to develop highly effective AMPs to decrease their working concentrations. However, the critical AMP concentration is most probably determined by the characteristic lengths of membrane deformations, which should weakly depend on the particular type of the membrane (bacterial or eukaryotic). In this view, the use of synergistic combinations of antimicrobial agents seems to be a more reliable alternative [24][25][26]. The synergism may be achieved by action of two peptides on two consequent stages of the bacteria damage. For example, the first peptide can permeabilize the plasma membrane thus allowing the second peptide to reach its target inside the cell [25]. In the work [26] the synergistic effect was shown to be based on that the binding of the first peptide (PGLa) to the membrane substantially facilitated the binding of the second peptide (magainin 2) via relaxation of unfavorable membrane curvature strain.
A special class of AMPs is so-called channel-formers, a typical member of which is gramicidin A (gA) [27,28]. This pentadecapeptide is produced by Bacillus brevis. In membrane environment gA attains β 6.3 -helical conformation [29,30]. The longitudinal axis of the helix is tilted but slightly with respect to the membrane normal; such an orientation is maintained by three tryptophan amino acids at gA C-terminal [28,31,32]. Two gA monomers located in the opposing leaflets of the membrane can reversibly form a cation-selective ion channel via transmembrane head-to-head dimerization [27,29]. Normally, the conducting dimer is stabilized by six hydrogen bonds established between amino acids at N-terminal of two gA monomers [33]. Both dimerization and dissociation of the dimer are suggested to pass through the common state of so-called coaxial pair (denoted as P), in which two gA monomers from opposing lipid monolayers are located one on top of the other [34][35][36]. The coaxial pair is considered to correspond to the top of the energy barrier of the dimerization/dissociation process, whereas the states of two monomers or conducting dimer correspond to (local) minima of the free energy. In all three characteristic states (far-separated monomers, coaxial pair, conducting dimer) gA induces elastic deformations of the membrane [31,32,[34][35][36]. In turn, elastic parameters of the membrane (thickness, spontaneous curvature, lateral tension, elastic rigidity) substantially affect the characteristics of the channel (lifetime and probability of formation) [37][38][39]. This allowed utilizing gramicidin as a molecular sensor of membrane elasticity [40][41][42]. Various membranedeforming inclusions, like amphipathic peptides, are also expected to influence the lifetime and probability of formation of the conducting channel via deformation-mediated lateral interactions with gA monomers, dimers and coaxial pairs. The interactions may differently contribute to the energies of these states, thus modifying the height of the energy barriers of gA dimerization/dissociation. Recently, we have demonstrated that even gramicidin monomers themselves may serve as such membrane-deforming inclusions, and they can influence the characteristics of the channel via deformation-mediated lateral interactions with the dimer and coaxial pair [36]. Lateral dimers of gA may be formed via connection of two gA monomers located in the same membrane leaflet by a linker [43][44][45]. In a series of works [43][44][45] it is demonstrated that the lifetime of so-called tandem channels of gA formed from two lateral dimers located in opposing membrane leaflets, exceed the lifetime of the regular gA channel by about three orders of magnitude. This effect can be also quantitatively explained by strong membrane-mediated lateral interactions between two transmembrane dimers of gA that comprise the tandem channel [35]. The experimentally observed strong influence of membrane-mediated lateral interactions on gA channel characteristics makes gA a convenient tool to validate predictions derived based on the theory of elasticity of lipid membranes.
The practical use of gA as an antimicrobial agent is very limited because of its high toxicity. In highly-organized organisms, relatively hydrophobic gA molecules incorporate into cellular plasma membranes without any significant selectivity towards bacterial cells [46,47]. Monomer-dimer equilibrium of gA is strongly shifted towards monomers. If gA concentrations are sufficient to raise the permeability of bacterial membranes, the permeability of eukaryotic membranes will raise as well. At smaller concentrations, the average number of gA channels is too small to cause any significant antimicrobial effect.
Here we attempt to build a rational basis for development of specific molecules amplifying the membrane-permeabilizing activity of AMPs in the target (bacterial) membrane. The amplifiers selectively binding to bacterial membranes are supposed to allow reducing working concentrations of AMPs and simultaneously increase their selectivity. Permeabilization of the membrane via formation of barrel-like pore, toroidal pore or ion channel is accompanied by elastic deformations of the membrane, and these deformations contribute to characteristics of the membrane-permeabilizing structures. In other words, the average conductivity of the membrane should depend on the spatial structure of its deformations. Rigid membrane inclusions that purposefully modify the deformation fields of the membrane can thus potentially regulate the activity of AMPs. The application of the amplifiers would allow separating the requirements to AMPs in terms of high specificity and activity. With the use of amplifiers, AMPs may have a relatively low selectivity, and the amplifiers may not have the ability to permeabilize the membrane at all. With the simultaneous addition of AMP in a low (non-toxic for eukaryotic cells) concentration and the amplifiers of its activity, AMP can be incorporated into membranes of almost all cells. However, the permeabilization should occur exclusively in bacterial membranes due to the selective incorporation of the amplifiers of AMP activity into bacterial membranes. Due to the selective incorporation into bacterial membranes of both the AMP itself and the amplifier of its activity, the effective distribution coefficient between the membranes of bacterial and eukaryotic cells of such a combined antimicrobial drug is determined by the product of the distribution coefficients of its individual components (AMP and activity amplifier).
In a number of works devoted to investigation of the molecular mechanisms of AMP action, a small number of peptides simultaneously acting on the membrane are considered [14,16]. However, the permeabilization of membranes by AMPs at a certain surface concentration is a stochastic process, and the average membrane conductivity induced by AMPs is a statistical value. Thus, it seems reasonable to use statistical approaches for understanding and adequate description of AMP action mechanisms and for estimation of the membrane conductivity under certain experimental conditions. For definiteness, here we develop the statistical model for conductivity of the membrane induced by gramicidin A, although the consideration can be generalized for other types of AMPs inducing elastic deformations. As potential amplifiers of the gA activity, we consider single-pass transmembrane peptides, e.g., WALP or KALP [48][49][50], and short amphipathic peptides. Both type of peptides (denoted as I, inclusion) are supposed to carry a chemical group responsible for the selective binding to bacterial membranes. We demonstrate that even low concentrations of amphipathic peptides (e.g., 1 peptide per 2500 lipids) and extremely low concentrations of short transmembrane peptides (e.g., 1 peptide per 100,000 lipids) are sufficient to substantially (several times) alter the average characteristics of the gA channel ensemble (i.e., lifetime and number of conducting dimers).

Algorithm of the Model Development
Consider the membrane with N 1 monomers of gA in the upper monolayer (denoted as U), N 2 monomers of gA in the lower monolayer (denoted as L), and m particles of the inclusion (denoted as I). At the moment, gA monomers from opposing membrane leaflets form k conducting dimers (denoted as D). The rate of transition of gA monomers between the leaflets is assumed to be very low, and thus the numbers of gA monomers in each leaflet N 1 , N 2 are constant. First, we derive the partition function Z 0 of the membrane with adsorbed gA molecules and incorporated inclusions. In the zeroth approximation, we consider the ideal system, where all possible membrane-mediated interactions of gA species and inclusions are absent. In order to account for the interactions of gA species and membrane inclusions (all particles with all ones), we utilize the Mayer cluster expansion [51][52][53]. In this approach, it is supposed that the configuration partition function of the non-ideal system can be expressed as a product of specific functions of pairwise interaction potentials; each function is averaged with respect to all possible positions of two interacting particles accounting for the Gibbs weight of the configurations. For gA monomers incorporated into two (upper and lower) leaflets of the membrane that are able to form conducting dimers, in the presence of inclusions, there are ten such pairwise interaction functions corresponding to the total number of ways to choose couples of particles from four types of particles present in the system: gA monomer in the upper monolayer (U), gA monomer in the lower monolayer (L), conducting dimer (D), inclusion I: U + U, U + L, L + L, U + D, L + D, U + I, L + I, D + D, I + I, D + I. In the limit of low concentrations, the pairwise functions are statistically independent and the average of their product can be substituted by the product of their average values. This substitution yields the first approximation of the partition function. Following van Kampen [52], we look for the higher-order corrections by means of correcting factors. In particular, the second correction should include three-particle correlations yielding the factors in the partition function in an amount equal to the number of ways to select three particles from four types of particles (U, L, D, I) present in the system; if we had only one kind of particles, there would be N × (N − 1) × (N − 2)/3! of such factors. Utilizing this approach, for known potential energy of particle pairwise interactions, one can calculate the coordination partition function in any desired order on concentrations. Generally, the factors for calculation of the partition function in the second order approximation on concentrations can be presented as the linear diagram (Figure 1a), third order-as the triangular diagrams (Figure 1b), fourth order-as square-like or fourth order diagrams shown in Figure 1c. In the diagrams, each vertex corresponds to the particle of the particular type, each edge corresponds to the function f AB , which characterizes the interaction of the particles A and B, standing in the corresponding vertices, f AB = exp[-w AB /(k B T)] − 1, where w AB is the pairwise interaction potential of particles A and B, k B is the Boltzmann constant, T is the absolute temperature. The interactions should be integrated over all possible configurations (positions) of the particles yielding the Mayer's cluster integrals, and then should be summed over all possible types of particles. All illustrated diagrams are irreducible [51][52][53]. Utilizing the described approach, for known potential energy of particle pairwise interactions, one can calculate the coordination partition function in any desired order on concentrations. In the present work, we consider the partition function in the fourth order on concentrations.

Pairwise Interaction Potentials
The profiles of pairwise interaction potentials calculated in the framework of the theory of elasticity of lipid membranes are illustrated in Figure 2 for configurations involving gA monomers. From Figure 2 it is seen that the interaction potentials with gA The potential energy of the pairwise interaction of the particles was considered as the energy of elastic deformations of the membrane in the corresponding configuration of the particles. The elastic energy was calculated in the framework of the theory of elasticity of lipid membranes [20,[34][35][36]54,55]. The gA monomer, dimer and membrane inclusion differed by the type of the boundary conditions they impose on the membrane deformations. The obtained potential energies allowed calculating numerically the Mayer's integrals of the necessary order. The Mayer's integrals were then used to obtain the partition functions of the system. The partition functions allowed calculating the dependences of the average lifetime of the gA channels as well as the average number of dimers in the membrane (i.e., the total membrane conductivity) on the concentrations of gA monomers and membrane inclusions (peripheral or transmembrane).

Pairwise Interaction Potentials
The profiles of pairwise interaction potentials calculated in the framework of the theory of elasticity of lipid membranes are illustrated in Figure 2 for configurations involving gA monomers. From Figure 2 it is seen that the interaction potentials with gA monomers are mostly repulsive except that for monomer-dimer (w UD ) and monomer-transmembrane inclusion (w UI ). The potentials decay rapidly with increasing distance between the particles, and at r ≈ 10 nm the interactions are negligibly weak.

Pairwise Interaction Potentials
The profiles of pairwise interaction potentials calculated in the framework of the theory of elasticity of lipid membranes are illustrated in Figure 2 for configurations involving gA monomers. From Figure 2 it is seen that the interaction potentials with gA monomers are mostly repulsive except that for monomer-dimer (wUD) and monomer-transmembrane inclusion (wUI). The potentials decay rapidly with increasing distance between the particles, and at r ≈ 10 nm the interactions are negligibly weak. Figure 2. Profiles of pairwise interaction potentials for configurations involving gA monomers. Red curve-wUU(r); green curve-wUL(r); blue curve-wUP(r); magenta curve-wUD(r); dark blue curve-wUI(r) (for transmembrane inclusion). The configurations are schematically illustrated in the right-hand-side of the figure. The membrane is shown as light gray stripe; solid black lines correspond to surfaces of the upper and lower monolayers; dotted line corresponds to the monolayer interface; gA monomers are shown as squares; the transmembrane inclusion is shown as the rectangle. The colors of the particles in the right-hand-side of the figure correspond to the colors of the interaction potential curves wAB. Red curve-w UU (r); green curve-w UL (r); blue curve-w UP (r); magenta curve-w UD (r); dark blue curve-w UI (r) (for transmembrane inclusion). The configurations are schematically illustrated in the right-hand-side of the figure. The membrane is shown as light gray stripe; solid black lines correspond to surfaces of the upper and lower monolayers; dotted line corresponds to the monolayer interface; gA monomers are shown as squares; the transmembrane inclusion is shown as the rectangle. The colors of the particles in the right-hand-side of the figure correspond to the colors of the interaction potential curves w AB .
All dependences w AB illustrated in Figure 2 have global minima. The minima are relatively deep when gA monomers interact with transmembrane particles: for dimer, w UD , and transmembrane inclusion, w UI , the depths of the minima are~3.5 k B T and~1.5 k B T, respectively ( Figure 2). The minima of the potentials that do not involve transmembrane particles are relatively shallow and their depths do not exceed about 0.5 k B T (curves w UU , w UL , w UP in Figure 2). In all configurations considered in Figure 2, the monolayer interface is not flat as the particle configurations and their induced deformations do not possess mirror symmetry with respect to the plane of the unperturbed membrane. Figure 3 illustrates the profiles of pairwise interaction potentials for configurations that do not involve gA monomers. From Figure 3 it is seen that the interaction potentials 7 of 26 in such configurations are mostly attractive. The potentials decay rapidly with increasing distance between the particles, and at r ≈ 10 nm the interactions are negligibly weak. All dependences w AB illustrated in Figure 3 have global minima, and their depths are quite large: from~4 k B T for w PI up to~11 k B T for w DD . All configurations illustrated in Figure 3 are mirror symmetric with respect to the plane of the unperturbed membrane; the monolayer interface is flat.
All dependences wAB illustrated in Figure 2 have global minima. The minima are relatively deep when gA monomers interact with transmembrane particles: for dimer, wUD, and transmembrane inclusion, wUI, the depths of the minima are ∼3.5 kBT and ∼1.5 kBT, respectively ( Figure 2). The minima of the potentials that do not involve transmembrane particles are relatively shallow and their depths do not exceed about 0.5 kBT (curves wUU, wUL, wUP in Figure 2). In all configurations considered in Figure 2, the monolayer interface is not flat as the particle configurations and their induced deformations do not possess mirror symmetry with respect to the plane of the unperturbed membrane. Figure 3 illustrates the profiles of pairwise interaction potentials for configurations that do not involve gA monomers. From Figure 3 it is seen that the interaction potentials in such configurations are mostly attractive. The potentials decay rapidly with increasing distance between the particles, and at r ≈ 10 nm the interactions are negligibly weak. All dependences wAB illustrated in Figure 3 have global minima, and their depths are quite large: from ∼4 kBT for wPI up to ∼11 kBT for wDD. All configurations illustrated in Figure 3 are mirror symmetric with respect to the plane of the unperturbed membrane; the monolayer interface is flat. Figure 3. Profiles of pairwise interaction potentials for configurations that do not involve single gA monomers. Red curve-wPI(r); green curve-wDI(r); blue curve-wPD(r); magenta curve-wDD(r); dark blue curve-wII(r) (for transmembrane inclusion). The configurations are schematically illustrated in the right-hand-side of the figure. The membrane is shown as light gray stripe; solid black lines correspond to surfaces of the upper and lower monolayers; dotted line corresponds to the monolayer interface; gA monomers are shown as squares; the transmembrane inclusions are shown as rectangles. The colors of the particles in the right-hand-side of the figure correspond to the colors of the interaction potential curves wAB.

Influence of Lateral Interaction on the Equilibrium Number of Dimers
The pairwise interaction potentials presented in Figures 2 and 3 were utilized to numerically calculate the Mayer's cluster integrals of the necessary order. The Mayer's integrals were then used to obtain the partition functions of the system. Some intermediate equations could not be solved analytically, and thus, to obtain the approximate solutions, we utilized the perturbation theory. The partition functions allowed obtaining the dependence of the dimer-monomer equilibrium constant on the concentration of gramicidin monomers, amphipathic membrane inclusions, and transmembrane inclusions in the membrane. This dependence leads to the shift of the equilibrium number of Figure 3. Profiles of pairwise interaction potentials for configurations that do not involve single gA monomers. Red curve-w PI (r); green curve-w DI (r); blue curve-w PD (r); magenta curve-w DD (r); dark blue curve-w II (r) (for transmembrane inclusion). The configurations are schematically illustrated in the right-hand-side of the figure. The membrane is shown as light gray stripe; solid black lines correspond to surfaces of the upper and lower monolayers; dotted line corresponds to the monolayer interface; gA monomers are shown as squares; the transmembrane inclusions are shown as rectangles. The colors of the particles in the right-hand-side of the figure correspond to the colors of the interaction potential curves w AB .

Influence of Lateral Interaction on the Equilibrium Number of Dimers
The pairwise interaction potentials presented in Figures 2 and 3 were utilized to numerically calculate the Mayer's cluster integrals of the necessary order. The Mayer's integrals were then used to obtain the partition functions of the system. Some intermediate equations could not be solved analytically, and thus, to obtain the approximate solutions, we utilized the perturbation theory. The partition functions allowed obtaining the dependence of the dimer-monomer equilibrium constant on the concentration of gramicidin monomers, amphipathic membrane inclusions, and transmembrane inclusions in the membrane. This dependence leads to the shift of the equilibrium number of conducting dimers, C D , as compared to the case of absent lateral interactions, C id D . Below, we graphically illustrate the dependence of the ratio C D /C id D = 1 + ∆C D /C id D (see Equation (24)) on the concentration of gramicidin monomers, amphipathic membrane inclusions, and transmembrane inclusions. For simplicity, we consider the case C U = C L = C 1 = C 2 corresponding to equally distributed gramicidin monomers over the membrane leaflets. Besides, we analyzed the consequences of membrane-mediated lateral interactions on gA analogues with N-terminal valine replaced by glycine ([Gly1]gA) or tyrosine ([Tyr1]gA). These analogues have significantly reduced the monomer-dimer equilibrium constants D 0 , with monomer-dimer equilibrium being shifted towards monomers as compared to gA [36]. The values of monomer-dimer equilibrium constants (product of dimerization constant and dimer lifetime) determined in the work [36] for gA, [Gly1]gA, [Tyr1]gA are as follows: 25,812 nm 2 , 250.8 nm 2 , 52.5 nm 2 , respectively.
The dependences of the equilibrium concentration of conducting dimers on concentrations of gramicidin monomers and membrane inclusions (amphipathic or transmembrane) are shown in Figure 4 for gA, [Gly1]gA, and [Tyr1]gA. The magenta curves in Figure 4 bound the range of concentrations where the correction to the main approximation calculated within the perturbation theory reaches 20%. This value of the correction was arbitrary chosen to qualitatively mark the limits of applicability of the perturbation theory.
Besides, we analyzed the consequences of membrane-mediated lateral interactions on gA analogues with N-terminal valine replaced by glycine ([Gly1]gA) or tyrosine ([Tyr1]gA). These analogues have significantly reduced the monomer-dimer equilibrium constants D0, with monomer-dimer equilibrium being shifted towards monomers as compared to gA [36]. The values of monomer-dimer equilibrium constants (product of dimerization constant and dimer lifetime) determined in the work [36] for gA, [Gly1]gA, [Tyr1]gA are as follows: 25,812 nm 2 , 250.8 nm 2 , 52.5 nm 2 , respectively.
The dependences of the equilibrium concentration of conducting dimers on concentrations of gramicidin monomers and membrane inclusions (amphipathic or transmembrane) are shown in Figure 4 for gA, [Gly1]gA, and [Tyr1]gA. The magenta curves in Figure 4 bound the range of concentrations where the correction to the main approximation calculated within the perturbation theory reaches 20%. This value of the correction was arbitrary chosen to qualitatively mark the limits of applicability of the perturbation theory.  From Figure 4 it is seen that increasing concentrations of gramicidin monomers (C U = C L ), amphipathic membrane inclusions and transmembrane inclusions lead to increase of the ratio C D /C id D . The most pronounced growth of the number of the conducting dimers is induced by transmembrane inclusions (Figure 4a-c): the doubling of the conducting dimer number is provided by only 2.5 particles per 1 µm 2 . Although the amphipathic membrane inclusions cause appreciable increase of the C D /C id D in substantially larger concentrations (e.g., hundreds of particles per 1 µm 2 , (Figure 4d-f)), practically, it is much easier to incorporate amphipathic peptides into the membrane as compared to transmembrane peptides.

Influence of Lateral Interaction on the Average Lifetime of Dimers
The membrane-deforming inclusions influence the lifetime of the conducting dimers of gramicidin. The dissociation constant, which is the inverse lifetime of the dimer, is determined by the difference of the energies of the coaxial pair and the dimer [34][35][36]. Partially, these energies include contributions from membrane elastic deformations, which, in turn, can be modulated by membrane-deforming inclusions. In order to characterize the effect of membrane-mediated lateral interactions, we consider the ratio τ/τ 0 of the average dimer lifetime, τ, and this value τ 0 in the "ideal" case of absence of any lateral interactions. The dependences of τ/τ 0 on concentrations of gramicidin monomers and membrane inclusions (amphipathic or transmembrane) are shown in Figure 5 for gA, [Gly1]gA, and [Tyr1]gA. The magenta curves in Figure 5 bound the range of concentrations where the interaction-related correction to the main "ideal" approximation calculated within the perturbation theory reaches 20%. This value of the correction was arbitrary chosen to qualitatively mark the limits of applicability of the perturbation theory. From Figure 5 it follows that increasing concentrations of amphipathic membrane inclusions and transmembrane inclusions lead to increase of the ratio τ/τ 0 . The most pronounced growth of the dimer lifetime is induced by transmembrane inclusions (Figure 5a-c): a 2-3 times increase of the lifetime is provided by only 2-3 particles per 1 µm 2 . Amphipathic membrane inclusions cause appreciable increase of the dimer lifetime in substantially larger concentrations-hundreds of particles per 1 µm 2 , (Figure 5d-f). The level curves in contour plots in Figure 5 are approximately parallel to the axis of gramicidin monomer concentration (C U ) meaning that the dependence of τ/τ 0 on C U is relatively weak. In the limit C U → 0 one can obtain the simple expression, Equation (31), for the dependence of the lifetime of the solitary dimer on the concentration of membrane inclusions. This dependence is illustrated in Figure 6 for transmembrane and amphipathic membrane inclusions.     From Figure 6 it is seen that the average lifetime of the solitary conducting dimer grows with increasing concentrations of transmembrane inclusions and amphipathic membrane inclusions. An about order-of-magnitude increase in τ/τ 0 is provided by 4.5 transmembrane peptides or 2400 amphipathic peptides per 1 µm 2 .

Discussion
In the present study we considered the consequences that the membrane-mediated lateral interactions have on characteristics of the gramicidin ensemble and solitary dimers. For the analysis, we developed the statistical model that accounted for the interactions in the fourth order on concentrations of membrane-incorporated particles. The interactions were assumed to be mediated by elastic deformations of the membrane, induced by the membrane-incorporated particles. The profiles of the potential energy of pairwise interactions of the particles were calculated in the framework of theory of elasticity of lipid membranes [20,[34][35][36]54,55]. Generally, the interaction involving gramicidin monomers are mostly repulsive (Figure 2), while the interactions involving transmembrane particles are mostly attractive (Figure 3). From the model it follows that gramicidin monomers themselves stimulate formation of conducting dimers excessive over the number of dimers would be in the "ideal" case of absent lateral interactions (Figure 4), while the average lifetime of the conducting dimer weakly depends on the monomer concentration ( Figure 5). Along with gramicidin species (monomers, dimers, coaxial pairs), we considered "inclusions"-membrane-deforming peptides, either transmembrane or amphipathic (peripheral). Short transmembrane peptides were predicted to substantially increase both the number of conducting dimers and the average lifetime of the dimers in very low concentrations, of the order of several peptide molecules per 1 µm 2 (Figures 4a-c, 5a-c and 6a). Amphipathic peptides alter the gramicidin characteristics less strongly: a substantial increase in the number of conducting dimers and the lifetime of the dimer requires concentrations of the order of hundreds-thousands of amphipathic peptides per 1 µm 2 (Figures 4d-f, 5d-f and 6b). Both transmembrane inclusions and amphipathic membrane inclusions act as amplifiers of gramicidin activity, as they stimulate formation of conducting dimers and increase the average lifetime of the conducting state (Figures 4-6).
If these inclusions have higher affinity towards particular membranes, e.g., they preferentially incorporate into negatively charged bacterial membranes, the inclusions would amplify the gramicidin activity predominantly in these membranes. The high selectivity of the inclusions and their amplification of gramicidin activity should allow reducing substantially the working concentrations of the gramicidin (e.g., below the toxic level for eukaryotic cells) still preserving a strong antimicrobial effect. Although the amphipathic membrane inclusions cause appreciable increase of the C D /C id D in substantially larger concentrations than the transmembrane inclusions do, practically, it is much easier to incorporate amphipathic peptides into the membrane as compared to transmembrane peptides. Amphipathic peptides themselves can cause formation of pores in the membrane. However, the critical concentrations of membrane disruption or lysis by amphipathic peptides such as magainin-2 are of the order of one peptide per 40-100 lipid molecules [1,8,10,13]. The typical concentration of 1000 amphipathic peptides per 1 µm 2 that is sufficient to provide substantial amplification of the gramicidin activity (Figures 4d-f and 6b) corresponds to approximately one amphipathic peptide per 2850 lipids, i.e., about two orders of magnitude smaller than the membrane lytic concentration.
The developed model is quite general; its possible applications are not limited to the case of gramicidin. The model may be utilized for analysis of amplifying or weakening effect of various membrane-deforming inclusions with respect to different types of AMPs. The only requirement is that at the sequential stages of the action, the AMP should induce different deformations of the membrane. For the example of gramicidin this means that the spatial distributions of elastic deformations of the membrane induced by gramicidin monomers, coaxial pair and dimer should be different. In the case of pore-forming AMPs, elastic deformations of the membrane can be induced by AMPs themselves, pores and membrane-incorporated inclusions. Accounting for the energy of the deformations, one can obtain the membrane-mediated interaction potentials, similar to those illustrated in Figures 2 and 3. The potentials can be further used to calculate the configuration partition functions and, consequently, to analyze the modulation effect of membrane inclusions on the pore-forming activity of AMPs, following the formalism developed in the "Materials and Methods" section.
The considered mechanism of the purposeful change in the equilibrium constant of the dimerization reaction is of a general nature and is not specific for gramicidin. In particular, we expect that in multiple systems, where membrane-deforming particles (monomers) can form higher-order aggregates (dimers, trimers, etc.) an addition of deformation-inducing inclusions can change the equilibrium distribution of the aggregate numbers. In such reactions, the inclusions do not play a role of a catalyst, as catalysts do not influence the chemical equilibrium. The inclusions strongly favor or stabilize a particular state or aggregate (e.g., the dimer in the case of gramicidin, Figure 3) thus shifting the equilibrium towards that state. This is rather similar to the action of the excessive external pressure on the equilibrium A + B ⇔ AB, where the volume of (A + B) reagents exceeds substantially the volume of AB product. In this case, a higher pressure favors smaller volume of the system, thus shifting the equilibrium to the right, towards AB product. The analogous effect should have an increased temperature if the direct reaction is endothermic, while the opposite reaction is exothermic. Commonly, in three-dimensional systems, the equilibrium constants weakly depend on the concentrations of the reagents. E.g., the Debye shielding works in solution; besides, when two oppositely charged particles combine, their charges cancel each other out and the combined particle (dimer) can interact with neighboring particles only via dipole and other interactions that are weaker than Coulomb interactions. In this sense, membranes are specific systems that allow particles to interact via relatively long-range deformations.
Fundamentally, the membrane with incorporated peptides can be considered as the system of many electric charges that mutually electrostatically interact. These interactions should alter when the membrane is disturbed by the inclusions, and thus the electrostaticbased interactions should contribute to the total energy of the system. However, when the membrane inclusion deforms the membrane, the bilayer thickness changes, the monolayer surfaces become curved, the monolayer surface can locally compress or stretch, and the axes of lipid molecules can deviate from the normal to the monolayer surface (i.e., the molecules can tilt). Let's focus on the bending (curving) of the monolayer. Denote the local curvature of the monolayer surface as J. Assuming J small, one can take Taylor series of the free energy F on the small deviation of the curvature from zero: Here, the first term is constant. The second term is zero, as we take Taylor series in the vicinity of the equilibrium state, where the free energy is minimal. In the third term, the is referred to as the bending modulus (denoted by B, see Materials and Methods section). In our model, we utilized the values of the elastic moduli that were determined experimentally. This means that our used value of the bending modulus includes (within the required accuracy) all possible contributions from all fundamental physical forces (electrostatic, van der Waals, entropy-based, etc.) that arise upon deviation of the monolayer surface curvature from zero. This is due to the fact that all these forces and all these contributions acted when the modulus was experimentally determined and contributed to its measured value. The same can be repeated for the other deformation modes (lateral stretching, tilt, etc.). Thus, in our calculations, the contribution from, e.g., electrostatic and van der Waals interactions acting in the system of many charges are accounted for within the required accuracy via the experimentally determined values of the elastic moduli.
For gramicidin A, the elastic approach for the description of monomer-dimer equilibrium was originally developed in 1986 [56] and then was modified in around 2000 to account for the effects of monolayer spontaneous curvature [57,58]. The general statement of the model is that the energy of elastic deformations arising in the vicinity of gA dimer depends quadratically on the hydrophobic mismatch (d 0 − l) between the length of the gA dimer (l) and bilayer thickness (d 0 ). Thus, the model is based on linear elastic theory (on the Hook's law), and is principally limited by the case of small deformations, i.e., (d 0 − l) should be much smaller than d 0 and l: 2(d 0 − l)/(d 0 + l) << 1. Note, that the particular limits of the linear theory applicability (i.e., the quantitative meaning of the sign "<<") cannot be determined within this theory; it can be only verified experimentally or, e.g., in molecular dynamics (MD). The fundamental flaw of the original hydrophobic mismatch-based elastic model was that it considered only the deformations induced by the gA dimer and assumed that neither gA monomer nor the gA coaxial pair (which we consider as the transition state in monomer-dimer reaction) deform the membrane. However, if the gA monomer length along the normal to the membrane surface is smaller than the monolayer thickness, elastic deformations should unavoidably appear in the gA monomer vicinity. As the gA monomer cannot impose boundary conditions onto the monolayer (bilayer) thickness, these deformations principally cannot be accounted for in the framework of hydrophobic mismatch-based models. This means, that the activation barriers of the dimerization and dissociation reactions of gramicidin are calculated not exactly correctly in such models, as the energy of deformations in the transition state and the energy of deformations induced by gA monomers are not taken into account. Thus, the single order parameter (bilayer thickness) introduced in the mismatch-based models is not sufficient to account for the deformations arising in the gA monomer vicinity. In our considerations, we utilize the additional order parameter, characterizing the average orientation of lipid molecules-socalled, director. The gA monomer was assumed to impose the boundary condition onto the director (Equation (38)) and this leaded to deformations of the membrane in the gA monomer vicinity.
Of note, the membrane-deforming inclusions alter the average energy of gramicidin monomers in the membrane. This may lead to shift of equilibrium distribution of monomers between the membrane and water phase. In this case, the actual concentrations of gramicidin monomers in the membrane should stand in the equations. Obviously, numerous fundamental physical forces act at the peptide-lipid boundary. In our model these interac-tion (electrostatic, van der Waals, entropy-related, etc.) determine the boundary conditions imposed on the parameters that describe the state of the lipid monolayer (director, thickness, etc.). The corresponding energy of interactions of peptide with the adjacent lipids yields a constant contribution. In our consideration we operate with interaction potentials that are normalized in a way to tend to zero as the distance between the peptides increases (Figures 2 and 3), and thus any constant contributions do not influence our model and its results.
It should be noted that the obtained dependences (Figures 4-6) are rather qualitative, since when calculating the interaction potentials, we used the boundary conditions derived from simple geometric considerations. The boundary conditions can be subjected to uncertainties, which can alter the quantitative values derived from the diagrams (Figures 4-6). However, our elastic model was successfully utilized in several works, and provided the results that were in quantitative (or, at least, semi-quantitative) agreement with the numerous experimental data and the data obtained in MD, including: the dependence of gramicidin A dimer dissociation probability on the lateral tension of the membrane; dependence of the probability of dimer formation on the lipid membrane thickness; dependence of the dimer lifetime on the spontaneous curvature of monolayers; the interaction energy profile of two gA dimers; the radial distribution of bilayer thickness around the gA dimer; the lifetime of the tandem gA channels; the characteristic time of flicker conductance of the tandem gA channels; the dependence of the characteristic time of gramicidin monomer-dimer equilibration on the gramicidin concentration, both normal dependence for gA, and anomalous dependence observed for gA derivatives [Gly1]gA and [Tyr1]gA [34][35][36]. Thus, the utilized parameters of the model and the model itself can be considered as reliable.
It is further worth noting that we used the large area limit, i.e., we assumed that although the concentrations of all particles are relatively low, the total number of the particles is still large that allow utilizing the statistical approach for the ensemble description. Practically, in the membrane of a relatively small area, the range of concentrations considered in our model can correspond to just a few particles. In this case, it is incorrect to apply the statistical description of the dimer-monomer equilibrium to obtain the average number of conducting dimers (Figure 4). However, it is still correct to use the developed model to describe the effect of membrane inclusions and gramicidin monomers on the lifetime of a solitary dimer ( Figure 6). Plasma membranes of both eukaryotic and bacterial cells may comprise a large number of membrane-deforming species (e.g., transmembrane and peripheral proteins, protein-lipid domains, etc.), which can take part in the modulation of the equilibrium between AMP monomers and cooperative membrane-permeabilizing structures. Thus, the introduced membrane-deforming inclusions would regulate the already modulated equilibrium. However, due to the high selectivity of the amplifiers towards the bacterial cells, they should nevertheless predominantly enhance the permeability of bacterial membranes.

General Expressions for Partition Functions
We develop the statistical model of gA ensemble in the presence of membranedeforming inclusions in the following way. Consider the membrane with N 1 monomers of gA in the upper monolayer (denoted as U), N 2 monomers of gA in the lower monolayer (denoted as L), and m particles of the inclusion. At the moment, gA monomers from opposing membrane leaflet form k conducting dimers (denoted as D). The rate of transition of gA monomers between the leaflets is assumed to be very low, and thus the numbers of gA monomers in each leaflet N 1 , N 2 are constant. The partition function of such system can be expressed as: where the summation is performed over the all possible numbers k of conducting dimers; the factor in the brackets is the number of ways to assemble k dimers from N 1 , N 2 monomers of gA; Z 0 (N 1 , N 2 , k, m) is the partition function of this system with fixed number of conducting dimers. In the zeroth approximation, we consider the ideal system, where all possible membrane-mediated interactions of gA species and inclusions are absent. In this approximation, the partition function can be written as follows: where is the Planck's constant; the temperature T is measured in thermal energy units; the first exponent is the Boltzmann factor (W D is the energy of single dimer); four factors in the brackets determine the phase space accessible by U, L, D, and I, respectively. These factors are proportional to TS/M or TS/M I , where S is the area of the membrane, M and M I are the effective masses of gA monomer or inclusion, respectively. The effective mass of the dimer is considered as a doubled effective mass of gA monomer; the effective mass is the actual mass of the peptide particle and the mass of lipid molecules involved in the peptide particle motion. In Equation (2) we omitted the internal partition functions of the particles (rotational, oscillatory, etc.) as they do not influence the results obtained. We introduce the function: where W is the energy of particular configuration of the system that depends on coordinates of all particles in the membrane; the integration is performed over all possible positions of the particles. In the case of non-interacting particles: (S is the area of the membrane). Utilizing the function Q, the partition function Equation (2) can be rewritten as: In order to describe the dependence of gA monomer-dimer equilibrium and the gA channel lifetime on gA monomer concentrations in the presence of inclusions, one should calculate the function Q taking into account the membrane-mediated lateral interactions of dimers, monomers, and inclusions. In order to account for the interactions of gA species and membrane inclusions (all particles with all ones), we utilize the Mayer cluster expansion [51][52][53] in up to the fourth order on particle concentrations, considering the interactions as pairwise. In this approach, it is supposed that the configuration partition function of the non-ideal system can be expressed as: where ϕ ij = exp(-w ij /T); w ij are pairwise interaction potentials of particles i and j; N is the total number of particles; the overline here and below denotes averaging with respect to all possible positions of the interacting particles. It can be explicitly demonstrated that ϕ 12 and ϕ 13 as well as ϕ 12 and ϕ 34 are statistically independent, and thus ϕ 12 ϕ 13 = ϕ 12 · ϕ 13 , ϕ 12 ϕ 34 = ϕ 12 · ϕ 34 . However, generally, ϕ 12 ϕ 13 ϕ 23 = ϕ 12 · ϕ 13 · ϕ 23 . The relation ϕ 12 ϕ 13 ϕ 23 = ϕ 12 · ϕ 13 · ϕ 23 can be used as the first approximation in the case of low concentration of the particles, as in this case the configurations when all three functions ϕ ij differ from unity are relatively rare. Applying this approach, we can express the function Q in the first approximation as follows: where ϕ AB = exp(-w AB /T) is the function corresponding to interaction of particles A and B. For gA monomers incorporated into two (upper and lower) leaflets of the membrane that are able to form conducting dimers, in the presence of membrane inclusions, there are ten such pairwise interaction functions corresponding to the total number of ways to choose couples of particles from four types of particles present in the system: gA monomer in the upper monolayer, gA monomer in the lower monolayer, conducting dimer, inclusion: U + U, U + L, L + L, U + D, L + D, U + I, L + I, D + D, I + I, D + I. In the limit of low concentrations, the pairwise functions are statistically independent and the average of their product can be substituted by the product of their average values. In the limiting case of large system (S → +∞, N → +∞, N/S = C = const) the expression Equation (7) can be rewritten as: where β 1,AB = (ϕ AB − 1)dr AB is the Mayer's first irreducible cluster integral [53], r AB = r A − r B is the vector directed from the particle B to the particle A. The power of the exponent in Equation (8) includes 10 terms, each of which corresponds to the certain pair of interacting particles. Following van Kampen [52], we look for the higher-order corrections by means of correcting factors. In particular, the second correction should include three-particle correlations yielding the factors in the partition function in an amount equal to the number of ways to select three particles from the four types of particles (U, L, D, I) present in the system; if we had only one kind of particles, there would be N × (N − 1)×(N − 2)/3! of such factors. Introducing the notation ϕ AB = 1+ f AB , within the requited accuracy one can obtain: where C ABC is the number of ways to choose particles of types A, B, or C from the total number of particle types present in the system; the summation is performed over all possible triplets of the particles. In particular, for gA monomers in the upper leaflet, C UUU = N U (N U − 1)(N U − 2)/3!; the number of ways to choose two gA monomers in the upper leaflet and the dimer is C UUD = N U (N U − 1)k/(2!·1!); the number of ways to choose gA monomer in the upper leaflet, gA monomer in the lower leaflet and the dimer is given by C ULD = N U N L k/(1!·1!·1!), etc.
In order to obtain the next correcting factor, one should consider all possible groups of four particles and calculate the factor ϕ 12 ϕ 13 ϕ 23 ϕ 14 ϕ 24 ϕ 34 D [52], where D is the approximation of the value ϕ 12 ϕ 13 ϕ 23 ϕ 14 ϕ 24 ϕ 34 obtained on the previous iteration. Generally, the factors for calculation of the partition function in the second order approximation on concentrations can be presented as the linear diagram (Figure 1a), third order-as the triangular diagrams (Figure 1b), fourth order-as fourth order diagrams shown in Figure 1c. In all diagrams, each vertex corresponds to the particle of the particular type, each edge corresponds to the function f. The interactions should be integrated over all possible configurations (positions) of the particles, and then should be summed over all possible types of particles (analogously to Equation (9)). All illustrated diagrams are irreducible [51][52][53]. Utilizing the described approach, for known potential energy of particle interaction, one can calculate the coordination partition function in any desired order on concentrations. We consider the partition function in the fourth order on concentrations. In this case, the expression for Q(N 1 , N 2 , k, m) can be presented as: where P is the polynomial (of the fourth order in our case) depending on the concentrations of particles of different types; the coefficients of P are determined by the diagrams illustrated in the Figure 1. Explicitly, the polynomial P can be written as follows:

Equilibrium between gA Monomers and Dimers in the Case of No Lateral Interactions
In the case when all possible membrane-mediated interactions of gA species and inclusions are absent, the partition function given by Equation (1) supplemented by Equations (2)-(5) can be calculated explicitly: where U is the Tricomi's (confluent hypergeometric) function. Formally, the obtained partition function Equation (12) allows solving the problem of determination of the average number of dimers. In particular, in the limit of the large area, it yields the well-known relation for the dimer concentration [27]: where C id D is the concentration of dimers in the ideal case of absent interactions; D 0 is the monomer-dimer equilibrium constant in this case. This relation directly follows from Equation (12): A similar relation can be analytically obtained directly from the partition function, Equations (1)-(5), by the saddle point method. Let's rewrite Equation (5) using the Stirling's formula: where const is some function independent on the number of dimers, k; f (k) is the logarithm of the expression standing under the summation sign in the upper line of the equation. The function f (k) has a maximum, which can be obtained from the equation: From this expression it is seen that in the limit of large number of particles, the maximum of f (k) is achieved at k = k 0 , where k 0 is determined by Equation (14). In the vicinity of the maximum: As N 1 , N 2 , k 0 are much larger than 1, the integral in Equation (15) is majorly accumulated in the vicinity of k 0 . This allows rewriting Equation (15) as: This partition function yields the same average number of dimers k = k 0 , as the partition function given by Equation (12). The square root in Equation (18) can be omitted, as the average number of dimers in the macroscopic limit is insensitive to this factor.

Equilibrium between gA Monomers and Dimers: The First Order of Interactions
Let's consider the partition function given by Equations (1)- (5) in the linear order on the number of dimers. In this case, the polynomial P(k), Equation (11), should be truncated on linear term on k: P = p 0 + α 1 k, where p 0 and α 1 are independent on k, although they can depend on N 1 , N 2 , m, T. Then the partition function can be expressed as (compare to Equation (12)): where W * D /T = W D /T − α 1 . Analogously to Equation (14), the average number k 1 of dimers in the first-order approximation can be obtained as: and, consequently, the concentration of dimers in the first-order approximation: where D 0 is the monomer-dimer equilibrium constant in the limit of low concentrations; the same D 0 stands in Equation (13).

Equilibrium between gA Monomers and Dimers: Higher-Order Terms
Below, we obtain the average number of dimers taking into account k 2 and higher order terms in Equation (11). Up to the fourth order, one can present the polynomial P(k) = p 0 + α 1 k + α 2 k 2 + α 3 k 3 + α 4 k 4 , where α 1 , α 2 , α 3 , α 4 are independent on k, although they can depend on N 1 , N 2 , m, T. In this case, Equation (16) yields: This equation cannot be solved analytically. To obtain the approximate solution, we utilize the perturbation theory. We assume that the average equilibrium concentration of dimers obtained taking into account k 2 and higher order terms in Equation (11), differ but slightly from the equilibrium concentration obtained in the first-order approximation as given by Equation (21), i.e., that for ∆C D = C D − C 1st D , the relation |∆C D | C 1st D holds. The Taylor series of Equation (22) up to the third order on C D yields: Further, we take Taylor series of C D with respect to C 1st D up to the fourth order, substitute C D into Equation (23) and under the assumptions C 1st D C 1 , C 2 express the coefficients of the Taylor series as: The series given by Equation (23) is formally correct when α 2 C dim0 1, α 3 C 2 dim0 1, α 4 C 3 dim0 1. Figures 2 and 3 demonstrate that all particles (U, L, D, I) on the average attract to transmembrane inclusions (gA dimers and transmembrane peptides), that results in positive values of corresponding Mayer cluster integrals β 1,AB , β 2,ABC , β 3,ABCE . The interaction energies are negative and their absolute values are quite large, of the order of several k B T. This allows assuming qualitatively, that the main contribution to the shift of monomerdimer equilibrium comes from interactions involving transmembrane inclusions, and thus all coefficients α i should be positive. Consequently, we obtain the lower estimate for the concentrations of dimers.

Lifetime of Dimers
Consider the gA dimer surrounded by gA monomers, other gA dimers, and membrane inclusions. The inverse average lifetime of the dimer, 1/τ, is determined by the product of the Boltzmann factor exp[-∆W/T] (where ∆W is the energy barrier of the dimer dissociation) and the frequency of attempts, ν, to dissociate: Following [34][35][36], we assume that the configuration of the coaxial pair of gA monomers corresponds to the top of the energy barrier. In this configuration, two gA monomers from opposing membrane leaflets stand one on top of the other and do not form the dimer. Thus, the energy barrier W is given by the difference of the energy of the pair and dimer configurations. The presence of gA monomers, gA dimers, and membrane inclusions influences the values of these two energies, and, consequently, alters the average lifetime of the dimer. In particular, in the works [43][44][45] it was demonstrated that the lifetime of tandem gA channels formed from two lateral dimers of gA, exceeds the lifetime of regular gA channels by about three orders of magnitude.
In order to obtain the inverse average lifetime of dimers in the membrane, one should average the expression Equation (25) with respect to positions of all particles present in the membrane: where W P,all is the energy of the system comprising gA monomers, dimers and inclusions supplemented by single coaxial pair of gA monomers; W D,all is the energy of the system comprising gA monomers, dimers and inclusions supplemented by single dimer (below, such dimer is denoted as D ). Actually, the relation (26) can be considered as a ratio of two partition functions. Utilizing the results on the formal structure of the partition functions obtained in the previous subsections, we can express 1/τ as follows: where τ 0 is the lifetime of single solitary dimer; K 1 , K 2 , K 3 are contributions from diagrams ( Figure 1) that include exactly one coaxial pair; H 1 , H 2 , H 3 are contributions from diagrams ( Figure 1) that include at least one dimer; ρ(k) = exp[f (k)] is the distribution function of the number of dimers. Below we present the explicit expressions for K 1 , H 1 , K 2 ; other contributions can be obtained in analogous way: Actually, it is not necessary to explicitly calculate the integrals standing in Equation (27), as in the macroscopic limit the corresponding integrals generally accumulate in the vicinity of the average value k = k , i.e., the distribution function ρ(k) is very narrow (close to the corresponding δ-function), and thus: In this equation, the correction to the dimer concentration calculated in accordance with Equation (24) should be taken into account. In the case of absence of membrane inclusions, one can explicitly check that in the first order on concentrations C 1 , C 2 from Equation (29) it follows: That coincides with the results obtained in the work [36]. Thus, in the regime of low concentrations, the main contribution to the lifetime correction is provided from interactions of monomers with pairs and dimers. Increasing concentrations lead to growing contributions from dimer-dimer and dimer-pair interactions. By considering the limit C 1 , C 2 → 0, one can obtain from Equation (29)

Calculation of Interaction Energies Based on Theory of Elasticity of Lipid Membranes
We calculate the energy of membrane-mediated interaction of gA monomers, dimers, pairs and membrane inclusions based on the theory of elasticity of lipid membranes. Each the particle induces elastic deformations in its vicinity. When the particles are far separated, their induced deformations are independent and their energies are additive. Upon mutual approach, the deformations overlap leading to effective lateral interaction between the particles. The method of calculation of the interaction energies is described in details in [34][35][36] for gA monomers and dimers, and in [20]-for amphipathic peptides. Here we briefly outline the algorithm of the calculation.
We introduce a Cartesian coordinate system Oxyz in such a way that the axis Oz is directed perpendicular to the plane of the undisturbed lipid bilayer; the coordinate origin O is placed at the monolayer interface. The upper monolayer lies in the half-space z > 0, the lower monolayer-z < 0. We describe elastic deformations of the lipid membrane considering it as a continuous elastic medium. The deformations are deemed small and their energy is calculated in quadratic order, analogously to the approach developed by Hamm and Kozlov [54]. The shapes of the upper and lower monolayer surfaces are characterized by z-coordinates of their points, described by functions H u (x, y) and H l (x, y), respectively. We describe z-coordinates of the monolayer interface by the function M(x, y). The average orientation of lipid molecules in the upper and lower monolayers is described by unit vectors n u and n l , respectively, called directors. Along with functions H u (x, y) and and r 0 = 1 nm, respectively [34][35][36]. The length and radius of the transmembrane inclusion were taken 2h I = 1.7 nm and r I = r 0 = 1 nm, respectively.

Numerical Minimization of the Elastic Energy Functional
The functional of the membrane elastic energy given by Equation (34) was minimized numerically, as the minimization cannot be performed analytically under the imposed boundary conditions. We utilized the finite element method with an adaptive mesh. The functions in the functional Equation (34) were replaced by their piecewise linear interpolants on each triangle of the computational mesh, the size of which decreased while approaching the membrane inclusions. The calculation was performed for meshes of different average sizes, and further the result was extrapolated towards the zero size mesh. The procedure of the numerical minimization of the elastic energy functional Equation (34) is described in details in the works [20,[34][35][36]65]. The resulting interaction potentials are presented in Figures 2 and 3.

Numerical Calculation of the Mayer's Cluster Integrals
The Mayer's cluster integrals β 1,AB , β 2,ABC , β 3,ABCE standing as coefficients in the polynomial, Equation (11), were calculated according to diagrams (Figure 1) utilizing the Monte Carlo method. For configurations of overlapping monomers in one monolayer, overlapping dimers, etc. the energy was set at the value equal to +∞. In the region where the interaction energies are finite, all calculated interaction energies were approximated by the function w(r) = e −qr P 8 (r), where q is some fitting constant, and P 8 (r) is the approximation polynomial of the eighth degree. The deviation of the approximation from the explicitly calculated energies did not exceed 2%. The explicit values of the integrals β 1,AB are presented in Table 1; β 2,ABC -in Table 2; β 3,ABCE -in Table 3.